%% Evaluating FOC for EB estimation (new)

function foc_val = foc(nu, hat_theta, rep_base_se2)
% INPUT: parameters nu, baseline study data hat_theta and rep_base_se2
% OUTPUT: evaluated value of FOC at given nu and tau_est

    [N, J] = size(hat_theta);
    
    % First part of the FOC
    first_term = 0;
    for i = 1:N
        for j = 1:J
            first_term = first_term + hat_theta(i,j)^2/((nu+rep_base_se2(i,j))^2) - 1/(nu+rep_base_se2(i,j));
        end
    end
    
    % Second part of the FOC
    second_term = 0;
    for i = 1:N
        % Calculate A_i, B_i, D_i, E_i
        A_i = 0;
        B_i = 0;
        D_i = 0;
        E_i = 0;
        
        for j = 1:J
            A_i = A_i + 1/(nu+rep_base_se2(i,j));
            B_i = B_i + hat_theta(i,j)/(nu+rep_base_se2(i,j));
            D_i = D_i + 1/((nu+rep_base_se2(i,j))^2);
            E_i = E_i + hat_theta(i,j)/((nu+rep_base_se2(i,j))^2);
        end
        
        second_term = second_term + D_i/A_i - 2*B_i*E_i/A_i + (B_i^2*D_i)/(A_i^2);
    end
    
    % Complete FOC value
    foc_val = first_term + second_term;
end